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1.  Summary 

In  this  document  we  report  on  activities  related  to  the  project  "New  physical  constraints  for  multi-frame  blind 
deconvolution",  grant  number  FA8655-12-1-2115,  in  the  second  year  of  funding.  Based  on  theoretical  work  on 
speckle  statistics  in  the  first  year  of  the  project  we  developed  a  deconvolution  algorithm  which  outperforms  state 
of  the  art,  while  being  free  of  parameters  which  normally  have  to  be  tweaked  by  the  user. 

In  what  follows,  we  summarize  all  relevant  developments,  also  some  from  Year  1,  and  give  examples  of  results 
obtained  on  the  AFRL's  3.5m  SOR  telescope  in  New  Mexico. 

2.  Introduction 

Images  of  space  objects  taken  with  ground-based  telescopes  are  blurred  by  atmospheric  turbulence.  The  problem 
can  be  alleviated  by  the  use  of  adaptive  optics  (AO)  and/or  image  post-processing.  Adaptive  optics  systems  work 
well  in  infrared  but  they  struggle  at  shorter  wavelengths.  With  AO  in  the  visible,  attainment  of  diffraction-limited 
resolution  is  currently  only  possible  in  conjunction  with  subsequent  deconvolution  (Figure  1).  An  optimal  way  to 
process  AO  images  is  the  one  based  on  the  knowledge  of  the  whole  imaging  process:  starting  from  the  influence 
of  turbulence,  through  the  optical  system  and  its  aberrations,  and  ending  with  (partial)  image  correction  by  AO. 
Such  a  physics-based  approach,  if  it  could  be  made  fully  automatic,  would  reduce  the  need  for  human 
intervention. 

Currently,  the  most  popular  approach  to  deconvolution  of  AO  images  is  the  so-called  "blind"  deconvolution  which 
relinquishes  the  afore-mentioned  determinism.  Instead,  it  is  assumed  that  smart  computer  algorithms  will  find 
the  optimal  way  to  deconvolve  the  data.  On  the  other  hand  blind  deconvolution  is  rarely  executed  blindly.  All 
available  methods  have  parameters  which  the  user  fine-tunes  until  the  most  visually-appealing  reconstruction  is 
achieved.  The  "art"  of  deconvolution  is  to  find  constraints  which  allow  for  the  best  estimate  of  an  object  to  be 
recovered,  but  in  practice  these  parameterized  constraints  often  reduce  deconvolution  to  the  struggle  of  trial  and 
error. 


* 


Fig.  1.  Images  of  Jupiter's  moon  Ganymede.  Left:  summation  of  thousand  short  exposures,  no  AO.  Middle: 
AO  switched  on.  Right:  AO  data  after  blind  deconvolution.  For  deconvolution  we  used  MISTRAL  package  [1]. 
Data  taken  with  the  AFRL's  3.5m  telescope  in  New  Mexico  at  Starfire  Optical  Range  (SOR),  wavelength  of  the 
observations  was  850nm,  AO  system  is  based  on  a  Shack-Hartmann  sensor  with  24  x  24  subapertures. 
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The  project  had  two  objectives:  (1)  to  develop  and  test  a  new  theory  of  speckle  statistics  based  on  the  model  of 
partially-developed  speckle,  and  (2)  to  test  if  inclusion  of  this  theory  in  deconvolution  removes  the  need  for 
human  intervention. 

3.  Methods,  Assumptions,  and  Procedures 

In  the  first  year  of  the  project  we  performed  a  large  set  of  experiments  involving  four  state-of-the-art 
deconvolution  algorithms.  The  setup  of  the  experiments  and  conclusions  were  reported  in  the  interim  report  and 
also  published  in  a  journal  paper  [2]  therefore  we  will  not  expand  on  them  here.  Worth  mentioning  is  that  the 
Adaptive  Wavelets  Maximum  Likelihood  Estimator  (AWMLE)  algorithm  developed  as  part  of  the  project  for  AO- 
assisted  observations  was  subsequently  re-written,  documented  and  packaged  for  public  release  in  2014. 

In  the  second  year  of  funding,  analytical  work  on  modelling  of  spatial,  temporal  and  spectral  integration  of 
speckle  statistics  was  carried  out.  Constraints  in  Fourier  domain  were  also  worked  out.  A  novel  methodology  of 
testing  the  reconstruction  fidelity  was  put  forward.  These  aspects  of  our  work  are  detailed  in  the  following 
sections. 

3.1  Integrated  speckle  statistics 

It  is  well  known  that  integration  (e.g.  adding  of  several  independent  realizations  in  one  pixel)  changes  speckle 
statistics  [3].  This  is  an  important  problem  in  imaging  through  turbulence  because  one  of  the  most  widely  spread 
solution  to  the  problem,  speckle  imaging,  relies  on  non-integrated  speckle  in  order  to  work.  On  the  other  hand, 
very  little  in  the  way  of  analytical  modelling  of  the  effect  of  integration  (for  example,  increasing  the  exposure  time 
of  the  observations)  has  been  attempted.  To  our  knowledge,  only  one  author  attacked  this  problem  but  the 
resulting  formulas  are  rather  complex  and  have  dependencies  on  several  parameters  which  cannot  be  obtained 
from  the  data,  like  phase  variance  [4].  We  aimed  to  provide  an  alternative,  simple  description  of  the  effect  of 
integration  on  various  speckle  statistics. 

Figure  2  illustrates  the  problem.  The  integration  parameter,  M,  is  small  in  the  top  row  and  large  in  the  bottom 
row.  M  can  be  thought  of  as  a  number  of  speckles  contributing  to  one  pixel.  While  the  change  in  M  due  to  spatial 
and  temporal  integration  has  been  the  subject  of  some  attention  in  the  context  of  scattering  of  rough  surfaces 
[3],  the  spectral  component  of  M  has  not,  and  this  constitutes  our  contribution.  Imagining  a  discrete  set  of 
increasing  wavelengths,  producing  radially  expanding  PSFs,  we  see  that  stepping  over  wavelengths  results  in  new 
speckles  traversing  a  pixel  under  observation  if  it  is  away  from  the  PSF  center.  The  further  a  pixel  is  from  the 
center  the  more  speckles  will  traverse  it  (compare  the  sharp  speckles  close  in  to  the  smooth  streaks  further  out 
from  the  center  in  the  rightmost  column  of  Figure  2). 

We  proposed  analytical  solutions  based  on  the  definition  of  M  (ratio  of  spatial,  temporal  or  spectral  auto¬ 
correlations  of  speckle  and  pixel  shape)  and  empirical  solutions  based  on  relation  of  M  to  first-order  speckle 
statistics: 


{h(x,y)) 

—7—7  »  =  fJVT 

var(h(x,yjj 
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(1) 


Distribution  A:  Approved  for  public  release;  distribution  is  unlimited. 


where  h(x,  y)  is  one  speckle  pattern,  Q  stands  for  mean  value  and  var(.)  stands  for  variance.  The  symbols  n,  v 

and  t  are  spatial,  spectral  and  temporal  integration  parameters.  The  first  two  need  only  be  estimated  once  for  a 
given  detector  and  filter.  This  is  the  approach  we  took  for  the  results  presented  here.  The  third  parameter  was 
estimated  by  computing  the  ratio  of  M  in  simulated,  noise-free,  instantaneous-exposure-time  simulations  to  M 
from  real  observations  at  SOR.  This  way,  it  was  found  to  be  around  5  for  our  datasets  [5].  An  illustration  of  proper 
choice  of  r  is  shown  in  Figure  3,  where  this  parameter  was  varied  around  5  to  test  the  effect  on  reconstruction  of 
one  particular  speckle  frame  from  SOR  observations.  This  can  be  considered  an  independent  test  for  our 
estimation  strategy  for  r.  The  methodology  pertaining  to  checking  the  fidelity  of  reconstructions  is  elucidated  in 
Section  3.3. 


Fig.  2.  Visual  effect  of  integration:  spatial  (left  column),  temporal  (middle  column)  and  spectral  (rightmost 
column).  For  each  case  a  red  square  was  used  to  denote  the  extent  of  one  pixel.  All  images  are  simulated.  In 
the  top  row  pixels  are  physically  small,  exposure  time  of  the  simulated  observations  is  very  short  (or  rather 
instantaneous)  and  the  band-pass  of  the  simulated  filter  is  negligible  (monochromatic  observations).  In  the 
bottom  row  the  situation  was  reversed:  the  pixel  size  is  large  compared  to  speckle  dimension,  as  many  as 
10  000  short  exposures  were  summed,  and  filter  band-pass  is  100  nm  (which  can  be  considered  large  for 
speckle  imaging). 


Fig.  3.  Reconstructions  of  one  particular  speckle  frame  (PSF)  based  on  various  choices  for  the  temporal 
integration  parameter  (spatial  and  spectral  integration  parameters  had  been  calculated  beforehand  and 
fixed).  The  rightmost  panel  shows  the  ground-truth  image.  Artificial  color  coding  and  square-root  scale. 
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3.2  Fourier-domain  constraint 


In  our  work  we  rely  on  rigorous  Bayesian  formulation  of  the  problem  of  finding  a  true  object  o  and  a  (stochastic) 
PSF  h  which  together  with  noise  generated  the  recorded  data  i: 


P(o,  h  1 1)  = - — -  (2) 

P(  i) 

The  theorem  states  that  the  conditional  probability  of  the  object  being  equal  to  reconstruction  o  and  the  PSF 
being  equal  to  reconstruction  h  given  that  we  recorded  data  i,  is  equal  to  the  product  of:  conditional  probability 
of  the  data  taking  on  the  value(s)  i  given  o  and  h,  probability  of  the  object  being  equal  to  o,  and  the  probability  of 
the  PSF  being  equal  to  h,  divided  by  the  probability  of  obtaining  the  data  i  which  is  always  taken  to  be  unity.  In 
the  maximum  a  posteriori  (MAP)  framework  one  finds  estimates  for  the  object,  6,  and  the  PSF(s),  h,  which 
jointly  maximize  p(o,h|i): 


[6,  h]  =  arg  max  p(o,  h|  i)  =  arg  max  /?(i|  o,  h)  x  p(o)  x  />(h)  (3) 

o,  h  o,  h 

Since  a  probability  density  function  p(o)  for  the  object  is  in  general  not  known,  solutions  in  the  form  of  functions 
with  desirable  mathematical  properties  (e.g.  noise  suppression,  edge  enhancement)  are  often  used.  As  these  ad 
hoc  formulas  are  not  real  probability  density  functions  (PDFs),  regularization  parameters  must  be  used  to  balance 
their  influence  on  the  cost  function  in  Equation  (3).  These  parameters  have  to  be  chosen  manually  [1].  Apart  from 
the  problem  of  choosing  the  right  value  for  them,  the  mere  form  of  the  prior  (e.g.  an  object's  power  spectral 
density)  could  be  applicable  only  to  a  limited  class  of  real  objects.  For  these  reasons  we  focus  on  PSF  statistics 
p(h).  This  part  of  Equation  (3)  was  almost  always  removed  from  the  MAP  approach  but  we  have  shown  how 
useful  it  is. 

In  Year  1  of  the  project  we  focused  on  the  gamma  distribution  for  the  focal-plane  intensity  statistics  of  h.  For  this 
purpose  we  developed  the  general  framework  for  inclusion  of  the  integration  parameter  for  intensity  (Section 
3.1).  We  also  proved  that  the  gamma  PDF  provides  a  good  model  for  the  integrated  speckle  (spatially,  temporally 
and  spectrally)  [6].  In  Year  2  we  also  developed  a  constraint  in  the  Fourier  domain  (where  Fourier  domain  means 
Fourier  transformed  focal-plane  intensity).  Under  the  assumption  of  normal  distribution  of  the  real  and  imaginary 
parts  of  the  optical  transfer  function  (OTF),  one  has  a  complete  model  for  their  joint  statistics: 
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_-M“) 

(9*)  =  e~  2 
(3)  =  0 , 


var(yi)  — 


2m(u) 


l  +  e-2D,(u)_2e-DM 


var(  3) 


1 


1  — e 


2m(u) 


2^(u) 


(4) 


where  D/uJ  is  the  phase  structure  function  evaluated  for  2-D  spatial  frequency  u  and  m(u )  is  the  number  of  OTF 
"cells"  at  that  frequency  (both  allowing  for  non-isotropic  turbulence)  [7]. 

This  prior  was  implemented  in  multi-frame  blind  deconvolution  algorithm.  Its  performance  is  illustrated  on 
simulated  observations  in  Figure  4  (real  l-band  speckle  PSFs  from  SOR  were  convolved  with  a  schematic 
representation  of  the  Hubble  Space  Telescope  which  acted  as  the  true  object). 


Fig.  4.  True  object  (leftmost  panel)  and  three  reconstructions  (from  left  to  right:  without  any  prior,  with  the 
spatial  constraint  gamma  prior,  and  with  the  Fourier  constraint  prior).  Refer  to  Figure  5  for  illustration  of  the 
blurred  input  data. 

3.3  Testing  reconstruction  fidelity 

Figures  4  and  5  are  typical  of  what  is  shown  in  most  image  processing  papers.  Reconstructions  of  the  objects  are 
successful  even  for  significant  blur  and  noise  corruption  (Figure  5,  bottom  row).  On  the  other  hand  we  are  putting 
constraints  on  the  PSFs,  not  on  the  objects,  and  therefore  it  is  more  interesting  to  check  the  accuracy  of  the  PSF 
reconstructions.  In  our  tests  we  fixed  the  object  and  we  tasked  the  algorithm  with  reconstructing  the  PSFs  only. 
This  way  we  could  check  the  benefit  of  applying  our  analytical  models  for  the  PSF  prior. 


We  allowed  the  algorithm  to  perform  20  000  iterations.  The  convergence  of  reconstructed  patterns  to  ground- 
truth  speckle  frames  was  clearly  visible  beyond  ca.  4000  iterations.  This  is  shown  in  Figure  6.  More  extensive 
testing  and  the  corresponding  results  are  discussed  in  Section  4. 
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Fig.  5.  From  left  to  right:  one  input  data  frame,  long  exposure,  and  reconstructions  obtained  with  our  MFBD 
code.  Top  row:  signal-to-noise  ratio  (SNR)  of  the  simulated  observations  was  100.  Bottom  row:  SNR  was  20. 
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Fig.  6.  Illustration  of  the  algorithm's  convergence.  The  code  was  supplied  with  perfectly  known,  fixed  object 
and  only  PSFs  were  reconstructed. 


4.  Results  and  Discussion 

In  Year  2  our  main  mode  of  checking  algorithm's  performance  was  by  looking  at  how  well  it  reconstructed  the 
ground-truth  PSFs.  We  were  interested  in  how  the  algorithm  performs,  given:  low-SNR  conditions,  bad  initial 
guess  for  the  average  PSF,  or  integration  parameter  M  in  the  gamma  model  for  the  intensity  prior. 

In  general,  it  was  observed  that  the  new  priors  significantly  helped  in  delaying  the  onset  of  noise  amplification. 
They  basically  provided  for  the  same  regularization  mechanism  as  traditional  object  priors,  with  the  major 
difference  that  PSF  priors  are  based  on  physical  models  of  the  aberrating  medium  (turbulence)  and  not  on 
properties  of  the  objects  being  imaged  (which  are  by  definition  unknown).  This  is  the  main  contribution  of  our 
work. 
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Figure  7  shows  how  a  lack  of  prior  in  traditional,  maximum-likelihood  deconvolution  leads  to  very  noisy 
reconstructions.  This  is  avoided  by  applying  a  gamma  prior  on  PSF  intensity  (central  panel). 
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Fig.  7.  PSF-prior  based  reconstruction  vs.  a  no-prior  execution.  In  both  cases  the  algorithm  was  supplied  with 
a  true  object  and  allowed  to  iterate  on  the  PSFs  (20  000  iterations). 

Another  important  aspect  of  our  work,  pursued  in  Year  1,  was  estimation  of  average  PSF  directly  from  the  target 
observations.  This,  in  principle,  allows  one  to  save  observing  time  (because  point-source  calibrators  are  not 
needed)  and  algorithm  execution  time  (because  the  trial-and-error  guessing  of  the  correct  first-guess  PSF  is 
eliminated).  In  order  to  remove  the  object  being  viewed  from  the  image  formation  equation  an  object  cancelling 
transformation  was  proposed  [8].  After  transforming  the  data,  information  about  the  object  cancels  out,  and  one 
is  left  only  with  moments  of  the  turbulent  OTF  which  is  customary  parameterized  by  the  Fried  parameter  r0.  In 
our  tests  this  method  achieves  accuracy  of  around  10%  on  r0  for  realistic  simulations  with  Gaussian  noise, 
broadband  filter  (AA/A  =  0.2)  and  only  100  input  frames  (Figure  8). 
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Fig.  8.  Demonstration  of  the  accuracy  of  the  "Fourier  contrast"  r0  estimation  method.  Based  on  simulated 
observations  of  artificial  satellites  with  a  0.5  m  telescope. 
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In  Year  2  we  checked  how  important  this  first  guess  for  the  average  PSF  is  for  successful  reconstructions  of  the 
individual  speckle  PSFs.  We  supplied  the  algorithm  with  r0  from  the  "Fourier  contrast"  method  [8,9]  which  was 
almost  equal  to  the  ground-truth  r0,  and  also  with  a  deliberately  mismatched  r0.  The  results  show  that  the  effect 
of  wrong  first  guess  becomes  especially  visible  at  low  SNR,  as  expected. 


Fig.  9.  Assuming  wrong  r0  has  little  impact  on  PSF  reconstruction  quality  when  SNR  of  the  dataset  is  high 
(here  100).  The  true  r0  was  10  cm. 


Fig.  10.  Incorrect  first  guess  for  r0  (and  correspondingly  for  the  average  PSF)  leads  to  badly  reconstructed 
PSFs  when  SNR  is  low  (here  20).  This  is  especially  visible  for  underestimated  r0  (leftmost  panel). 


5.  Conclusions 

In  our  project  "New  physical  constraints  for  multi-frame  blind  deconvolution",  grant  number  FA8655-12-1-2115, 
we  have  developed  new  models  for  statistics  of  speckle  images.  The  models  are  based  on  a  theory  which  is 
fundamentally  different  from  the  approaches  developed  in  the  1970s  for  speckle  imaging.  We  make  use  of  the 
equations  developed  for  scattering  off  rough  surfaces  ("partially-developed  speckle").  The  models  were  checked 
against  the  baseline  theory. 
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This  theory  was  employed  in  the  task  of  restoring  the  quality  of  images  taken  through  turbulence.  In  blind  tests, 
the  new  multi-frame  blind  deconvolution  provided  significantly  better  image  quality  than  one  of  the  codes  used 
by  AFRL  personnel  at  SOR.  Our  MFBD  code  can  be  executed  autonomously,  i.e.  without  tweaking  of  parameters 
and  without  specifying  a  first-guess  PSF.  The  project  has  resulted  so  far  in  one  journal  paper  (two  more  in 
preparation)  and  five  papers  in  conference  proceedings. 
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7.  List  of  Symbols,  Abbreviations,  and  Acronyms 

AO  -  adaptive  optics 
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AWMLE  -  adaptive  wavelets  maximum-likelihood  estimator 


MFBD  -  multi-frame  blind  deconvolution 
PDF  -  probability  density  function 
PSF  -  point-spread  function 
SNR  -  signal-to-noise  ratio 
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